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1.  Summary 

Under  the  auspices  of  a  one  year  contract  with  the  Air  Force  Research  Laboratory,  research  was 
conducted  on  advanced  imaging  of  underground  structures  using  artificial  electromagnetic 
sources  such  as  the  High  frequency  Active  Auroral  Research  Program  (HAARP)  Ionospheric 
Research  Instrument.  Non-invasive  imaging  of  underground  structure  is  important  for  the 
detection  of  hidden  tuimels  and  other  hazards,  as  well  as  resource  exploration,  mineral 
exploration,  and  environmental  contamination  problems.  Under  this  contract,  we  explored  the 
following  areas  related  to  underground  tunnel  detection:  1)  imaging  algorithms  including 
minimum  structure  inversions  and  parametric  inversions,  2)  analysis  of  the  sensitivity  and 
resolution  of  surface  electromagnetic  data  to  subsurface  tunnels,  3)  analysis  of  data  collected 
over  known  tunnel  structures,  both  by  the  US  Geological  Survey  and  by  APTI,  Inc.,  and  4) 
modification  of  our  robust  processing  algorithms  to  deal  with  AMT  data  and  for  controlled 
source  data. 

While  our  minimum  structure  imaging  codes  had  previously  been  written  and  tested,  we 
developed  new  imaging  codes  that  performed  parametric  inversions,  which  are  inversions  for 
simple  geometric  models  rather  than  fully  inhomogenous  models.  Two  such  codes  were 
developed,  one  that  inverted  data  for  sharp  boundaries,  and  one  that  inverted  data  for  a 
rectangular  tunnel  in  a  homogeneous  background  media.  These  inversion  algorithms  are  useful 
for  testing  the  possible  classes  of  models  consistent  with  a  given  set  of  data.  The  imaging 
algorithms,  along  with  simple  forward  modeling,  were  used  to  examine  the  sensitivity  and 
resolution  of  surface  EM  data  to  subsurface  tunnels.  Our  analysis  indicated  that  tuimels  needed  to 
be  at  approximately  a  depth  to  diameter  ratio  of  3:1  jDr  less  to  be  resolvable.  Actual  field  data 
collected  over  known  tuimels  showed  that  in  a  realistic  situation,  however,  that  even  this  may  be 
difficult  to  obtain  because  of  low  AMT  signal  levels,  and  geologic  noise.  A  phased  array 
transmitter  similar  to  HAARP  was  used  in  one  test  across  the  Silver  Fox  mine  in  Alaska,  but  the 
usefulness  of  the  transmitter  in  this  situation  was  limited  at  best.  Higher  power  transmitters  will 
likely  be  more  useful.  Finally,  our  robust  processing  algorithms  have  been  adapted  to  deal  with 
some  of  the  problems  with  AMT  data  and  controlled  sources,  and  they  show  considerable 
promise  for  future  work. 
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2.  Introduction 

Non-invasive  imaging  of  natural  and  man-made  underground  structure  continues  to  be  an 
important  issue  in  natural  resource  exploration  (oil,  gas,  minerals,  and  geothermal), 
environmental  contamination,  and  underground  hazards  detection  (tunnels  and  unexploded 
ordnance,  for  example).  Sub-surface  imaging  can  be  accomplished  using  a  variety  of  geophysical 
data  such  as  seismic,  gravity,  magnetic,  and  electrical  measurements.  In  particular,  electrical 
geophysical  techniques  (e.g.,  d.c.  resistivity,  magnetotellurics,  and  electromagnetic  (EM)  . 
soundings)  are  valuable  in  subsurface  exploration  because  the  bulk  electrical  conductivity  of  the 
subsurface  is  diagnostic  of  parameters  such  as  the  amount  and  connectivity  of  pore  fluids,  the 
presence  of  voids,  the  presence  of  chemical  contaminants,  and  the  presence  of  anomalously 
conductive  fluids.  Even  when  they  do  not  directly  detect  the  presence  of  fluids,  contaminants,  or 
voids,  electrical  methods  can  provide  valuable  information  on  structural  parameters  and  soil 
properties.  For  example,  an  image  of  the  sub-surface  electrical  conductivity  can  be  related  to 
geologic  and  structural  features  that  would  indicate  oil  and  gas  deposits  and  therefore  be  useful 
information  for  resource  exploration. 

Perhaps  one  of  the  most  common  and  widely  used  of  electrical  geophysical  methods  is 
the  exploration  method  known  as  magnetotellurics  (MT)  (Vozoff,  1991).  In  the  MT  method, 
measurements  of  the  naturally-occurring  electric  and  magnetic  fields  at  earth's  surface  are  used  to 
image  the  subsurface  electrical  conductivity.  The  frequency  range  spaimed  by  MT  measurements 
is  typically  from  500-0.001  Hz,  but  sometimes  higher  frequency  measurements  are  made  up  to 
20,000  Hz  in  the  method  known  as  audio-frequency  MT,  or  AMT.  Recent  advances  in  hardware 
and  MT  data  processing,  modeling,  and  interpretation  have  led  to  many  successes  in  crustal 
geologic  and  geophysical  exploration.  Nonetheless,  MT  still  suffers  from  certain  inherent 
limitations  including  low  signal  strength  in  the  1.0  to  0.1  Hz  region  (‘dead-band’)  and  around  1 
to  2  kHz  (the  AMT  dead-band),  variable  signal  levels,  and  man-made  electrical  noise  that 
obscures  the  natural  field  variations  and  leads  to  poor  data  quality. 

To  overcome  these  limitations,  geophysicists  have  employed  artificial  sources  to  generate 
large  amplitude  EM  fields  that  can  be  measured  in  a  maimer  similar  to  the  traditional  MT 
method.  This  technique,  a  variant  of  the  MT  method,  is  known  as  controlled-source  MT 
(CSMT),  and  the  sources  are  typically  long  grounded  dipoles  or  large  loops  of  wires  that  are 
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energized  with  alternating  currents  to  induce  EM  fields  in  the  earth.  CSMT  methods  overcome 
the  problems  with  low  signal  levels,  noise,  and  the  dead-hand,  but  unless  very  large  sources  and 
field  generators  are  used,  the  penetration  depths  for  CSMT  are  limited  to  a  few  kilometers  at 
best.  Additionally,  computer  simulations  of  EM  fields  in  a  heterogeneous  earth  for  finite  sources 
such  as  groimded  dipoles  is  time-consuming  and  computationally  intensive,  whereas  modeling 
EM  fields  from  a  plane  wave  source,  such  as  in  the  MT  method,  is  much  quicker  and  less 
computationally-intensive.  Even  though  geophysicists  often  go  to  great  lengths  in  CSMT  surveys 
to  approximate  plane  wave  fields  (e.g.,  making  measurements  far  away  fi'om  the  source),  there 
are  always  source  effects  in  the  measurements  that  can  make  CSMT  data  difficult  to  interpret. 

Therefore,  a  controllable  EM  somce  that  could  generate  signals  that  appeared  locally  as 
plane  waves  and  whose  signals  were  significantly  stronger  than  background  levels  would  be  of 
special  interest  for  subsurface  imaging  since  it  would  overcome  the  limitations  of  the  traditional 
MT  and  CSMT  methods.  The  phased  array  radio  transmitter  of  the  High  frequency  Active 
Auroral  Research  Program  (HAARP)  has  the  potential  to  be  just  such  a  controlled  plane-wave 
source.  By  controlling  the  modulation  frequency  of  the  transmitter,  which  was  originally 
designed  for  upper  atmospheric  and  solar-terrestrial  research,  EM  waves  can  be  generated  at  the 
frequencies  needed  for  geophysical  exploration  of  the  shallow  subsurface.  However,  the 
transmitted  power  levels  that  are  presently  sustainable  at  the  HAARP  facility,  which  is  only  in 
the  development  phase  of  its  construction,  are  much  lower  than  the  power  levels  planned  when 
future  phases  are  complete.  Consequently,  the  transmitted  signals  are  only  observable  within  the 
near  vicinity  of  the  facility,  but  when  at  full  power,  it  is  possible  that  the  transmitted  signals  will 
be  sufficiently  strong  to  be  used  for  exploration  and  imaging  at  greater  distances. 

The  potential  advantage  of  using  a  source  like  the  HAARP  transmitter  is  that  since  the 
observed  signals  would  appear  locally  as  plane  waves,  one  could  use  the  much  simpler  and  faster 
analysis  and  modeling  codes  developed  for  the  MT  method  and  not  have  to  worry  about  the 
effects  of  finite-dimension  sources.  Additionally,  EM  waves  generated  by  the  HAARP 
transmitters  would  not  suffer  from  the  ‘dead-hand’  limitations  of  the  MT  natural  fields,  and  in 
fact  may  allow  the  collection  of  interpretable  data  even  in  areas  with  high  levels  of  man-made 
electrical  noise.  Finally,  since  the  HAARP  transmitters  have  the  potential  to  generate  EM  waves 
at  low  frequencies,  one  can  explore  deeper  crustal  structure  that  are  not  accessible  to  the  usual 
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CSMT  sources.  The  viability  of  using  the  HAARP  transmitters  as  a  controlled-source 
geophysical  exploration  tool  has  already  been  demonstrated  (MacEnany,  1 997,  personal 
communication)  in  tests  carried  out  by  Advanced  Power  Technologies,  Inc.  (APTI)  for 
underground  tunnel  detection. 

In  this  report,  we  summarize  research  carried  out  on  the  applicability  of  using  HAARP 
for  undergroimd  tunnel  detection.  First,  we  outline  the  general  concepts  behind  our  imaging 
algorithms  and  extension  to  parametric  models.  Second,  we  discuss  the  sensitivity  of  EM 
measurements  to  undergroimd  tunnels  using  simple  synthetic  models  and  two-dimensional 
modeling  and  inversion  algorithms.  Third,  we  present  an  analysis  of  data  collected  by  the  both 
the  US  Geological  Survey  and  APTI  over  known  tuimel  structures.  Finally,  we  summarize  our 
recent  advances  with  robust  processing  codes. 

3.  Inversion  Algorithms 

The  goal  of  EM  imaging  (also  known  as  tomography  or  inversion)  is  to  infer  the  electrical 
conductivity  of  the  earth's  subsurface  on  the  basis  of  measurements  of  the  electric  and  magnetic 
fields  at  the  surface  or  in  boreholes.  The  problem  can  be  expressed  as  an  inverse  problem  of  the 
form 

d  =  A(jn)  +  e . 

The  data  vector  d  has  the  apparent  resistivities  and  phases  for  the  TM  and/or  TE  modes  for  all 
frequencies  and  locations.  The  resistivity  of  the  earth  is  parameterized  by  a  model  vector  m,  A  is 
a  functional  that  predicts  the  theoretical  value  of  d  for  a  particular  m,  and  e  is  the  observational 
error  in  d.  In  general,  the  earth's  resistivity  (p)  varies  significantly  both  vertically  and  laterally. 
Therefore,  it  is  appropriate  to  consider  w  to  be  a  function  of  3-D  position.  However,  owing  to 
the  computational  demands  of  3-D  model  simulations  and  because  2-D  approximations  are  often 
adequate,  in  this  report  we  will  consider  m  to  be  a  function  of  2-D  position.  To  enforce  positivity 
(p>0)  it  is  convenient  to  let  m  be  the  logarithm  of  p.  For  numerical  implementation,  the  model 
function  m  is  discretized  in  terms  of  a  2-D  grid  of  rectangular  blocks. 

The  functional  A  is  defined  implicitly  by  Maxwell’s  equations.  At  MT  frequencies, 
conduction  currents  dominate  over  displacement  currents.  Thus,  the  physical  process  is  one  of 
EM  induction,  with  EM  fields  diffusing  through  the  conductive  earth  rather  than  propagating  as 
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EM  waves.  Except  for  simple  classes  of  resistivity  models,  it  is  necessary  to  solve  Maxwell’s 
equations  numerically.  Numerical  techniques  applied  to  the  MT  problem  include  Madden’s 
(1972)  transmission  network  method,  the  finite  difference  method  (Mackie  et  al.,  1994),  and 
finite  element  method  (Waimamaker  et  al,  1984),  among  others. 

The  main  problem  with  electromagnetic  tomography  is  that  it  is  a  very  ill-posed  inverse 
problem.  Electromagnetic  measurements  can  be  made  only  at  a  finite  number  of  locations 
(surface  or  borehole),  whereas  the  conductivity  of  subsurface  materials  is  a  function  of  position 
which  can  vary  strongly  in  all  three  dimensions  and  therefore  has  an  infinite  number  of  degrees 
of  freedom.  Moreover,  the  inverse  problem  is  strongly  nonlinear  and  it  is  not  possible  to  “back- 
project”  or  migrate  the  data  directly  into  a  meaningful  image  of  the  subsurface,  as  is  done  with 
seismic  reflection  data.  An  iterative  process  is  therefore  necessary.  Finally,  the  fact  that  the  EM 
waves  diffuse  through  the  media  rather  than  propagate  as  waves  limits  the  basic  information 
content  of  the  data. 

3.1  Minimum  structure  inversion 

Because  of  the  ill-posedness  of  the  EM  inverse  problem,  it  must  be  stabilized  by  applying 
regularization  techniques.  Our  2-D  MT  inversion  algorithm  is  based  on  the  method  of  Tikhonov 
regularization  (Tikhonov  and  Arsenin,  1977).  This  method  defines  a  ‘regularized  solution’  of  the 
inverse  problem  as  a  model  that  minimizes  an  objective  flmctional  of  the  form 

d>  =  (4/  -  A(jn)y  !  +T  S{m) . 

Here,  t  is  a  positive  number  known  as  the  ‘regularization  parameter’.  The  first  term  of  O  is  a 
weighted  error  sum  of  squares  (‘misfit’)  for  the  model  m.  For  each  datum,  a  is  the  standard 
deviation  of  the  error  e.  The  second  term  is  t  times  a  functional,  S(m),  known  as  the  ‘stabilizing 
functional’. 

There  is  no  consensus  in  the  geophysical  community  about  the  best  choice  of  the 
stabilizing  functional,  S.  Some  workers,  motivated  by  a  stochastic  inversion  framework  (e.g. 
Tarantola,  1987),  define  S  in  terms  of  a  prior  model,  and  prior  covariance  operator,  R^m- 
This  approach  has  been  taken  by  Park  and  Van  (1991)  and  Zhang  et  al  (1995)  in  the  3D  d.c. 
resistivity  problem.  Alternatively,  many  workers  choose  5  as  a  measure  of  the  spatial  roughness 
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of  m,  which  is  typically  associated  with  the  size  of  its  first  or  second  spatial  derivatives.  The 
motivation  in  this  case  is  to  find  the  simplest  explanation  of  the  data  in  the  form  of  models  with 
minimum  structure  (Constable  et  al,  1987;  Smith  and  Booker,  1988;  Rodi,  1989).  The  minimum 
structure  approach  has  been  used  in  3D  resistivity  by  Ellis  and  Oldenburg  (1994),  while  several 
workers  have  employed  it  in  electromagnetic  inversion  problems  (Jiracek  et  al,  1987;  deGroot- 
Hedlin  and  Constable,  1990;  Smith  and  Booker,  1991;  Newman,  1995).  However,  the  precise 
definition  of  5  differs  among  the  workers  who  have  taken  this  approach. 

Our  current  tomography  algorithm  provides  different  options  for  the  stabilizing 
functional,  but  in  general  we  have  found  a  very  simple  functional  involving  only  second 
derivatives  of  m  to  be  effective: 

S{m)=  \(y^mfdx 

where  denotes  the  Laplacian  operator.  We  point  out  that  one  carmot  take  for  granted  that  the 
minimization  of  O  is  well-posed  for  any  choice  of  stabilizing  functional.  Unless  S  constrains 
derivatives  of  m  of  sufficiently  high  order,  numerical  solutions  of  the  inverse  problem  will 
depend  strongly  on  the  model  discretization  (grid)  and  the  model  will  contain  grid  artifacts 
regardless  of  how  fine  the  grid  is.  Numerical  experiments  we  have  performed  show  evidence  of 
this. 

Another  reason  for  choosing  minimum  structure  algorithms  is  related  to  the  physical 
process  that  controls  EM  fields  in  conductive  media.  At  the  frequencies  of  interest  for  geologic 
exploration,  EM  fields  diffuse  through  conductive  media  rather  than  propagate  as  waves.  An 
analogy  is  heat  conduction  through  a  metal  plate.  Consequently,  it  is  difficult  to  resolve  sharp 
boimdaries  in  the  absence  of  other  information.  Although  models  with  more  complicated 
structure  and  sharp  boundaries  may  be  consistent  with  the  observed  data,  they  are  rarely  required 
by  the  data,  and  this  distinction  is  an  important  issue  in  determining  what  features  are  resolvable 
by  the  data. 

3.2  Optimization  Algorithms 

EM  tomography  is  a  computationally  intensive  problem.  The  forward  problem  entails  solving  a 
large  matrix  system,  and  numerical  optimization  algorithms  are  needed  to  perform  the 
minimization  of  the  objective  functional,  <D. 
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Our  current  tomography  software  includes  two  algorithms  for  minimization  of  the 
objective  functional  (Rodi  and  Mackie,  1999).  One  is  the  “successive  linearized  inversion” 
scheme  described  by  Zhang  et  al  (1995)  and  based  on  earlier  work  by  Mackie  and  Madden 
(1993).  This  scheme,  which  can  be  categorized  as  a  version  of  the  Gauss-Newton  method 
(Tarantola,  1987),  solves  a  linearized  inverse  problem  in  an  iterative  loop  until  an  acceptable 
value  of  O  is  achieved.  We  can  then  solve  the  linearized  problem  iteratively  using  a  conjugate 
gradient  (CG)  technique.  The  CG  technique  requires  the  computation  of  linearized  forward  and 
adjoint  operations;  Zhang  et  al  developed  highly  efficient  algorithms  to  perform  these  operations, 
building  on  earlier  work  by  Madden  and  Mackie  (1989),  Rodi  (1976)  and  Mackie  and  Madden 
(1993).  The  EM  forward  problem  is  itself  solved  with  another  version  of  the  CG  method. 

The  other  optimization  algorithm  we  have  implemented  is  a  nonlinear  conjugate  gradient 
algorithm  (Rodi  and  Mackie,  1999),  an  approach  also  used  by  Ellis  and  Oldenburg  (1994). 

Unlike  the  Gauss-Newton  method,  nonlinear  CG  does  not  solve  a  linearized  inverse  problem  at 
each  step  of  an  iteration.  Instead,  at  each  step  it  performs  a  line  (one-parameter)  minimization 
along  a  given  direction  in  model  space.  The  model  direction  at  each  step  is  derived  firom  the 
gradient  of  the  objective  functional  for  the  current  model  and  the  direction  and  gradient  from  the 
previous  step.  A  pre-conditioning  operator  is  applied  to  the  gradient  vector  in  an  attempt  to  find 
the  most  important  directions  in  model  space  in  the  early  iterations.  The  same  efficient 
algorithms  for  forward  and  adjoint  operations,  used  in  the  other  inversion  scheme,  may  also  be 
used  in  nonlinear  CG. 

Our  experiments  with  the  nonlinear  CG  method  to  date  indicate  that  it  is  extremely  fast  if 
a  good  pre-conditioning  operator  is  used.  For  example,  inversions  of  real  data  with  up  to  120 
stations  and  20  frequencies  per  station  were  accomplished  on  a  desktop  workstation  in  less  than 
an  hour  with  the  nonlinear  CG  scheme  whereas  the  same  inversion  would  not  have  been  possible 
with  standard  matrix  inversion  Gauss-Newton  algorithms.  The  nonlinear  CG  method  has  been 
tested  on  a  wide  variety  of  real  and  synthetic  data  and  has  been  found  to  give  exemplary  results. 

3.3  Parametric  Inversions 

Minimum  structure  inversions  are  useful  because  they  yield  models  that  contain  only  the  minimal 
amount  of  structure  in  the  subsurface  required  to  fit  the  observed  data,  and  they  do  not  contain 
features  that  are  incidental  to  the  way  in  which  the  models  were  obtained.  Additionally,  it  is 
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possible  to  test  various  geological  h)rpotheses  in  the  minimum  structure  framework  by  the 
inclusion  of  an  a  priori  model.  The  inversion  routine  then  finds  the  minimum  variations  away 
from  the  a  priori  model  that  are  necessary  to  fit  the  data. 

However,  there  may  be  situations  where  smooth  gradational  models  are  not  optimal  in 
terms  of  geological  and  geophysical  interpretation.  For  example,  there  are  many  geological 
provinces  where  the  subsurface  geology  is  fairly  simple  and  well  known — it  may  consist  of 
basically  two  or  three  distinctive  units  which  are  separated  by  sharp  boundaries.  Another 
example  is  the  case  of  a  turmel  in  an  otherwise  fairly  homogeneous  background.  Here,  the  tunnel, 
which  is  very  resistive,  is  sharply  distinct  from  the  more  conductive  background.  The  problem, 
then,  is  to  determine  the  geometry  of  those  units  in  a  particular  location  within  that  geological 
province.  With  minimum  structure  inversions,  one  would  only  be  able  to  determine  gradational 
representations  of  those  boundaries  rather  than  the  boundaries  themselves.  While  this  may  be 
consistent  with  the  physics  of  EM  diffusion  in  a  conductive  media,  it  is  not  ideal  for  the  person 
who  needs  to  know  the  geometry  of  these  boundaries. 

One  way  to  deal  with  this  is  by  use  of  parametric  inversion  algorithms  (Portniaguine  and 
Zhdanov,  1995).  Instead  of  inverting  for  the  smoothest  model  possible,  one  inverts  for  only  a  few 
parameters  relating  to  simple  geometrical  bodies.  For  example,  one  could  specify  to  invert  for  the 
parameters  of  a  circular-shaped  body  in  the  subsurface  that  best  fit  the  observed  data.  The 
inversion  would  then  find  the  best-fitting  location,  radius,  and  conductivity  contrast  for  this  body. 
Although  such  structures  would  not  necessarily  be  required  by  the  data,  as  in  a  minimum 
structure  inversion,  it  would  allow  one  to  test  how  well  simple  sub-surface  geometries  fit  the 
observed  data. 

In  our  research,  we  have  developed  two  parametric  inversion  algorithms:  one  algorithm 
that  inverts  for  sharp  boundaries  with  variable  geometry,  and  another  that  simply  inverts  for  a 
rectangle  in  a  homogeneous  background  (the  tunnel  problem).  In  both  algorithms,  it  is  necessary 
to  derive  a  scheme  that  projects  discrete  boundaries  and  resistivities  onto  the  finite  difference 
mesh  used  for  the  forward  modeling.  For  the  sharp  boundary  inversion,  we  assign  discrete 
boundaries  by  specifying  nodes  with  straight  lines  connecting  the  nodes.  The  resistivity  of  the 
layers  are  allowed  to  vary  by  assigning  resistivity  nodes  within  the  layer  and  assuming  the 
resistivity  varies  linearly  between  nodes.  When  interfaces  cut  across  the  blocks  in  the  finite 
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difference  mesh,  we  use  a  simple  scheme  that  assigns  the  resistivity  of  that  block  as  a  weighted 
average  of  the  resistivity  above  and  below  the  interface.  Sensitivities  to  changes  in  both  the 
interface  node  locations  and  the  resistivities  of  the  layers  can  be  computed  efficiently  by  the 
calculus  chain  rule  and  use  of  the  same  efficient  forward  and  adjoint  operations  involved  in  the 
minimum  structure  inversion  scheme.  For  the  tunnel  inversion,  we  assume  a  homogeneous 
background  and  simply  project  the  tunnel  onto  the  finite  difference  mesh,  again  using  a  simple 
weighted  average  of  resistivities  where  the  tuimel  boundary  crosses  blocks  in  the  finite  difference 
mesh. 

3.4  Simple  Example  of  Parametric  versus  Smooth  Inversion 

Synthetic  MT  data  were  generated  for  a  simple  two-layer  model  with  a  sloping  ramp  interface, 
but  the  details  of  the  model  were  unknown  until  after  imaging  with  the  sharp  boundary  inversion 
algorithm  (the  model  and  synthetic  data  were  kindly  provided  by  a  colleague).  Random  Gaussian 
noise  at  a  level  of  5%  were  added  to  the  synthetic  data  and  these  data  were  then  imaged  by  both 
the  minimum  structure  algorithm  (Figure  1)  and  the  sharp  boundary  algorithm  (Figure  2).  It  was 
later  determined  that  the  sharp  boundary  inversion  came  within  a  few  percent  of  the  correct 
location  and  geometry  of  the  interface.  It  is  clear  from  the  figures  that  although  the  minimum 
structure  inversion  results  in  a  model  that  has  the  correct  overall  features,  it  would  be  difficult 
from  that  result  to  locate  boundaries  between  the  conductive  and  resistive  imits.  In  this  simple 
case,  the  sharp  boimdary  inversion  is  superior  to  the  smooth  model  result.  We  should  stress 
again,  however,  that  the  sharp  boundaries  are  not  inherently  resolvable  from  the  data  except 
when  we  specify  a  priori  that  the  resulting  model  will  have  sharp  boundaries. 

An  example  using  the  second  parametric  inversion  routine,  to  solve  for  a  simple 
rectangular  tunnel  in  a  homogeneous  backgroimd  media,  is  illustrated  in  the  next  section. 
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Smooth  model  inversion  of  Synthetic  Ramp  Data 
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Figure  1.  Smooth  model  inversion  of  synthetic  ramp  data. 
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4.  Sensitivity  of  EM  data  to  underground  tunnels 


4.1  Geometric  and  Resistivity  Effects 

Of  primary  importance  in  detecting  underground  tunnels  is  the  sensitivity  of  surface  EM  data  to 
the  size  and  depth  of  xmderground  structures.  That  is,  tunnels  that  are  too  deep  or  too  small  will 
not  have  an  observable  surface  signature  and  therefore  will  not  be  detectable  regardless  of  the 
imaging  algorithm  used.  Consequently,  it  is  critical  to  the  success  of  underground  tunnel 
detection  to  determine  the  depths  and  dimensions  at  which  underground  tunnels  have  some 
reasonable  chance  of  exhibiting  a  surface  signature,  above  and  beyond  the  normal  background 
geologic  signature. 


Figure  3.  TM  mode  amplitude  response  of  2x2  meter  tunnel  at  2  meters  depth.  The 
apparent  resistivity  responses  are  shown  for  various  resistivities  of  the  tunnel. 
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Therefore,  we  have  examined  the  sensitivity  of  EM  data  to  subsurface  tunnels  by 
computing  synthetic  data  for  simple  models,  using  finite-difference  modeling  algorithms.  In 
particular,  we  have  concentrated  on  models  that  attempt  to  simulate  subsurface  turmels  at 
different  depths  and  sizes.  The  first  model  tested  was  a  2m  by  2m  square  (of  varying  resistivity) 
buried  in  a  homogeneous  10  ohm-m  backgroxmd  media.  This  is  perhaps  the  simplest  model  and 
represents  an  upper  boundary  for  detection  of  subsurface  tunnels  (the  real  earth  is  not 
homogeneous  and  field  measurements  always  are  contaminated  by  noise).  Figure  3  shows  the 
TM  mode  apparent  resistivity  as  a  function  of  frequency  and  tunnel  resistivity  for  this  model 
when  the  top  of  the  simulated  tunnel  is  at  a  depth  of  2  m.  This  result  is  for  a  station  0. 125  m  off 
the  center  of  the  simulated  tunnel.  Figure  4  shows  the  TM  mode  phase  response  for  the  same 
model  and  also  for  vaiying  tunnel  resistivities.  Over  most  of  the  frequency  range,  the  TM 
response  appears  as  a  static  shift. 


Figure  4.  TM  mode  phase  response  of  2x2  meter  tunnel  at  2  meters  depth.  The  phase 
responses  are  shown  for  various  resistivities  of  the  tunnel. 
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The  TE  mode  apparent  resistivity  and  phase  for  this  same  model  are  shown  in  Figures  5 
and  6.  These  figures  demonstrate  a  well-known  aspect  of  MT,  which  is  that  TE  mode  data  are 
sensitive  to  conductive  bodies  and  have  very  little  sensitivity  to  resistive  bodies.  This  is  due  to 
the  gathering  of  current  into  the  infinite  extent  (in  2D  models)  of  a  conductive  body,  whereas  the 
currents  merely  flow  around  a  resistive  body.  From  these  figures  one  can  deduce  that  a 
conductive  body  is  more  resolvable  than  a  resistive  body,  and  indeed  this  is  borne  out  by 
inversions  of  these  synthetic  data.  Figure  7  shows  the  minimum  structure  inversion  result  for  the 
synthetic  data  where  the  subsurface  body  was  1  ohm-m,  whereas  Figure  8  shows  the  result  where 
the  body  was  1 000  ohm-m.  The  true  model  is  a  2m  by  2m  square  at  a  depth  extent  from  2  to  4m. 
Note  the  conductive  body  is  well-resolved,  because  of  the  strong  TE  mode  response,  whereas  the 
resistive  body  is  less  well-resolved,  because  it  only  has  weak  TM  mode  response  and  no  TE 
mode  response. 


Figure  5.  TE  mode  amplitude  response  of  2x2  meter  tunnel  at  2  meters  depth.  The 
apparent  resistivity  responses  are  shown  for  various  resistivities  of  the  tunnel. 
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This  demonstrates  that  even  under  ideal  conditions,  tuimels,  because  they  are  likely  to  be 
resistive,  will  be  difficult  to  detect. 

The  synthetic  data  shown  in  the  previous  examples  were  computed  using  our  own  finite 
difference  modeling  algorithms.  To  ensure  against  numerical  inaccuracies,  we  have  compared 
our  finite  difference  responses  to  the  responses  computed  for  the  same  models  using  Phil 
Wannamaker’s  finite  element  code.  The  comparison  is  shown  only  for  the  cases  of  the  1000 
ohm-m  square  and  1  ohm-m  square  in  a  10  ohm-m  background.  The  results,  Figures  9  and  10, 
indicate  excellent  agreement  between  the  two  modeling  algorithms,  with  less  than  0.5% 
differences. 


Figure  6.  TE  mode  phase  response  of  2x2  meter  tunnel  at  2  meters  depth.  The  phase 
responses  are  shown  for  various  resistivities  of  the  tunnel. 
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Conductive  model  inversion  results 
TM  and  TE  mode 


Figure  7.  The  minimum  structure  inversion  of  synthetic  data  over  a  conductive  (I  ohm-m) 
2x2  meter  tunnel  at  a  depth  of  2  meters  in  a  background  media  of  10  ohm-m.  Note  the  body 
is  well  resolved. 
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Resistive  model  inversic 

TWI  and  TE  mode 


Figure  8.  The  minimum  structure  inversion  of  synthetic  data  over  a  resistive  (1000  ohm-m) 
2x2  meter  tunnel  at  a  depth  of  2  meters  in  a  background  media  of  10  ohm-m. 
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Apparent  Resistivity  (ohm-m) 


Figure  9.  TM  mode  amplitude  response  comparing  our  finite  difference  algorithm  and  Phil 
Wannamaker’s  finite  element  algorithm. 
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TM  mode  Response  of  2m  tunnel,  2m  in  depth 
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Figure  10.  TM  mode  phase  response  comparing  our  finite  difference  algorithm  and  Phil 
Wannamaker’s  finite  element  algorithm. 


Next,  we  varied  the  depth  to  the  top  of  the  2m  by  2m  simulated  tunnel,  and  in  each  case, 
computed  the  synthetic  responses  for  varying  resistivities  of  the  tunnel.  The  TM  mode  apparent 
resistivity  is  shown  for  the  cases  when  the  depth  to  the  top  of  the  tmmel  is  4m  (Figure  1 1)  and 
6m  (Figure  12).  The  TE  mode  apparent  resistivity  for  the  case  when  the  depth  to  the  top  of  the 
tunnel  is  6m  is  shown  in  Figure  13.  It  is  clear  from  these  figures  that  by  the  time  the  subsurface 
body  is  at  a  depth/diameter  ratio  of  approximately  3-4,  the  response  is  extremely  weak  and  would 
be  very  difficult  to  detect  under  normal  field  conditions.  The  other  important  point  from  these 
simulations  is  that  the  frequencies  that  are  sensitive  to  observable  subsurface  tunnels  are  typically 
above  1000  Hz,  a  range  where  it  is  sometimes  extremely  difficult  to  get  good  data  because  of  the 
AMT  dead  band,  an  issue  that  will  be  discussed  shortly. 
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Figure  11.  The  TM  mode  amplitude  response  of  a  2x2  meter  tunnel  at  4  meters  depth.  The 
apparent  resistivities  are  shown  for  different  tunnel  resistivities. 
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Figure  12.  The  TM  mode  amplitude  response  of  a  2x2  meter  tunnel  at  6  meters  depth.  The 
apparent  resistivities  are  shown  for  different  tunnel  resistivities. 
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Figure  13.  The  TE  mode  amplitude  response  of  a  2x2  meter  tunnel  at  6  meters  depth.  The 
apparent  resistivities  are  shown  for  different  tunnel  resistivities. 


The  previous  examples  showed  the  results  of  using  minimum  structure  inversion 
algorithms  to  invert  the  synthetic  data  from  simple  models  simulating  shallow  turmels.  For  the 
case  where  the  tunnel  had  a  depth/diameter  ratio  of  1 ,  the  minimum  structure  algorithms  were 
able  to  resolve  both  conductive  and  resistive  bodies  quite  accurately.  We  then  applied  our 
parametric  inversion  algorithms  to  these  same  data,  solving  for  the  best  location  of  a  rectangular 
body  as  well  as  the  resistivity  of  the  body  and  the  background.  The  model  resulting  from  a 
parametric  inversion  of  synthetic  data  for  the  conductive  body  is  shown  in  Figure  14,  and  the 
model  resulting  from  a  parametric  inversion  of  synthetic  data  for  the  resistive  body  is  shown  in 
Figure  15.  In  both  cases,  the  results  from  the  parametric  inversion  match  nearly  perfectly  the 
parameters  for  the  true  synthetic  model  (2mx2m  v^ith  the  top  of  the  tunnel  at  2m  depth). 
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Conductive  model  inversion  results 

JM  and  TE  mode 


Figure  14.  Parametric  inversion  for  a  conductive  tunnel  of  dimensions  2x2  meters  at  2 
meters  depth. 
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Resistive  model  inversion  results 
TM  and  TE  mode 


resistivity  (ohm-m) 


Figure  15.  Parametric  inversion  for  a  resistive  tunnel  of  dimensions  2x2  meters  at  2  meters 
depth. 
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4.2  The  AMT  Dead  Band 


Natural  electromagnetic  fields  arise  from  a  variety  of  sources.  At  low  frequencies,  the  fields  arise 
mainly  from  the  interaction  of  earth’s  magnetic  field  with  the  solar  wind.  At  frequencies  between 
10  and  500  Hz,  the  fields  are  due  primarily  to  lightning  strikes  that  produce  EM  waves  that  travel 
around  in  the  earth-ionosphere  waveguide.  In  this  spectrum  of  natural  fields,  there  are  two  well 
known  dead  bands  (1  to  0.1  Hz,  and  1000  to  2000  Hz),  that  arise  from  a  lack  of  signal  and  from 
attenuation  in  the  earth-ionosphere  cavity.  The  most  severe  of  these  dead  bands  is  the  AMT  dead 
band  in  the  1 000  to  2000  Hz  range.  This  dead  band  is  important  for  tunnel  detection  because  it  is 
right  in  the  range  that  is  usually  critical  for  detecting  subsurface  tunnels.  We  have  carried  out 
independent  research  on  this  topic  and  found  that  even  using  24  bit  digitizing  systems,  in  general 
the  magnetic  field  spectrum  in  the  AMT  dead  band  is  below  the  inherent  system  noise  of 
commercial  coil  magnetometers.  Additionally,  we  have  found  that  the  natural  fields  in  this  band 
tend  to  be  of  greatest  magnitude  during  the  night  hours.  What  this  means  in  practice  is  that 
magnetic  field  measurements  in  this  frequency  band  are  very  difficult  to  make  reliably  and  hence, 
great  care  must  be  taken  to  ensure  accurate  measmements.  This  difficulty  is  obvious  in  some  of 
the  data  analyzed  later  in  this  report.  Electric  field  data  do  not  suffer  this  same  problem,  because 
the  electric  fields  can  now  be  measured,  using  24  bit  systems,  without  any  preliminary 
electronics  or  filtering.  The  main  impact  of  these  findings  is  that  a  reliable  source  would  be 
extremely  useful  for  collecting  AMT  data  for  tuimel  detection. 

5.  Analysis  of  Data 

In  this  section,  we  analyze  four  sets  of  MT  data  that  were  collected  over  known  tunnel  structures. 
The  first  three  data  sets  were  collected  by  Doug  Klein  from  the  USGS  and  were  provided  in 
order  to  test  competing  inversion  algorithms  for  accuracy  and  robustness.  The  fourth  data  set  was 
collected  by  APTI  over  the  Silver  Fox  mine  in  Alaska. 

5. 1  San  Xavier  Data 

The  data  in  this  experiment  were  collected  along  a  profile  oriented  338  degrees  relative  to  true 
north.  Electric  and  magnetic  data  were  collected  every  10  meters  along  this  profile.  An  EMI  ‘TX- 
IMl’  transmitter  located  approximately  300  m  southwest  of  the  line  was  used  as  a  controlled 
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inductive  source.  The  purpose  for  using  a  controlled  source  in  geophysical  surveys  is  to  obtain 
data  with  high  signal/noise  ratios,  especially  in  environments  where  natural  MT  signals  may  be 
difficult  to  collect  because  of  local  electrical  noise  sources. 

The  supplied  data  files  from  the  blind  test  were  input  into  a  WinGlink'  data  base  for 
further  analysis  and  consisted  of  both  controlled-source  data,  and  some  natural  field  MT  data  that 
were  collected  at  only  a  few  of  the  sites.  A  visual  assessment  of  the  data  indicated  generally  poor 
data  quality  with  a  wide  scatter  between  frequencies.  An  example  of  a  typical  sounding  is 
presented  in  Figure  16,  which  shows  the  sounding  curves  from  site  BA17. 


Figure  16.  An  example  of  a  typical  sounding  curve,  site  BA17. 


Before  interpretation,  the  data  were  manually  edited  on  a  site-by-site  basis,  first  eliminating 
obvious  outliers,  then  using  the  D+  algorithm  (which  computes  the  most  consistent  apparent 


'  WinGlink  is  a  data  base  and  interpretation  program  for  magnetotelluric  and  other  geophysical  methods  commonly 
used  in  geophysical  exploration,  and  was  developed  by  Geosystem,  srl. 
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resistivity  and  phase  cur\'es  for  a  given  set  of  data)  to  help  eliminate  inconsistent  data  points. 
While  this  is  a  partially  subjective  exercise,  it  does  result  in  data  that  are  less  noisy. 

After  editing  the  data  from  each  site,  the  data  were  rotated  into  alignment  with  the  profile 
direction  and  then  output  for  inversion  using  our  minimum  structure  inversion  algorithm.  We 
started  the  inversions  from  a  uniform  halfspace  of  30  ohm-m.  A  minimum  structure  inversion  of 
these  data  is  shown  in  Figure  17.  We  must  stress  that  since  the  data  are  of  generally  poor  quality 
and  only  marginally  useful,  conclusions  based  on  the  results  of  the  inversion  should  therefore  be 
taken  with  appropriate  caution.  The  txmnel  under  this  line  is  located  at  approximately  site  BAIO. 
The  inversion  puts  in  a  conductive  feature  at  this  location,  which  according  to  Doug  Klein  is  due 
to  fluid-filled  fractures.  We  would  normally  expect  a  tunnel  to  exhibit  a  resistive  signature  rather 
than  a  conductive  signature.  However,  it  is  possible  that  the  inversions  are  inaccurate  because  of 
the  poor  data  quality. 


Figure  17.  Minimum  structure  inversion  of  San  Xavier  MT  data. 


5.2  SSC  tunnel 


The  data  for  this  experiment  were  collected  at  the  SSC  tunnel  in  Waxahachie,  Texas  by 
Doug  Klein  of  the  U.S.  Geological  Survey.  AMT  data  were  collected  at  15  stations  spaced 
approximately  10  meters  apart  along  a  profile  line  perpendicular  to  the  timnel.  The  timnel  itself  is 
believed  to  have  been  approximately  5  meters  in  diameter  at  a  depth  of  50  meters,  giving  a 
depth/diameter  ratio  of  10.  The  AMT  data  were  collected  using  the  MT-1  system  manufactured 
by  EMI.  Natural  source  AMT  data  were  collected  at  each  site;  additionally,  lower  frequency  MT 
data  were  collected  at  a  few  of  the  sites,  and  controlled-source  data  were  collected  at  3  sites. 
Overall  data  quality  is  poor,  similar  to  that  observed  in  the  previous  data  set. 

All  acquisition  runs  at  each  site  were  combined  into  one  site  response,  converted  to  an 
EDI  file,  and  then  input  into  a  WinGlink  data  base.  The  responses  at  each  site  were  rotated  to  the 
strike  direction  of  162  degrees  magnetic,  and  the  TM  (transverse  magnetic)  and  TE  (transverse 
electric)  modes  were  output  for  inversion.  Inversions  were  carried  out  starting  from  a  uniform  10 
ohm-m  halfspace.  Several  inversions  were  nm,  but  the  main  features  were  consistent  among  the 
various  runs.  An  example  of  one  such  run  is  shown  in  Figure  18.  A  resistive  feature  is  imaged, 
but  at  a  depth  shallower  than  the  expected  tunnel.  This  is  likely  an  artifact  of  the  poor  data 
quality. 
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Figure  18-  Minimum  structure  inversion  of  SSC  data. 

5.3  uses  Blind  Experiment  #2 

The  data  from  the  USGS  blind  experiment  #2,  supplied  by  Doug  Klein  of  the  USGS, 
were  supplied  as  impedance  values,  and  these  were  converted  to  apparent  resistivity  and  phases 
and  then  used  for  generating  smooth  2D  models.  The  data  were  synthetic  data,  generated  for  a 
model  with  3  tunnels  at  depth  to  diameter  ratios  of  3:1, 10:1,  and  greater.  We  first  inverted  the 
noise-free  synthetic  data  using  our  minimum  structure  inversion  algorithm,  with  the  results 
shown  in  Figure  19.  There  is  an  indication  of  a  shallow  resistive  feature  located  near  the  4th 
station  (just  to  the  right  of  the  50m  mark  on  the  Figure)  and  at  a  depth  of  around  15-20  m.,  w'^hich 
is  likely  the  shallow-  tunnel.  The  tunnels  located  at  greater  depths  appear  not  to  be  resolvable 
from  this  data  set. 
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Inversion  of  USGS  Blind  Dataset  2 
TM  and  TE  data,  fit  to  1.0  sigma 
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Figure  19.  Minimum  structure  inversion  of  noise-free  USGS  data  from  Blind  Experiment 

#2. 


To  simulate  geologic  noise,  we  then  added  5%  random  Gaussian  noise  to  the  synthetic  data,  and 
then  imaged  the  data.  The  result,  shown  in  Figure  20,  is  similar  to  the  previous  result  in  that  only 
1  tunnel  is  identified,  but  is  probably  more  realistic  than  inverting  noise-ffee  data. 
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Inversion  of  USGS  Blind  Dataset  2 


TM  and  TE  data,  fit  to  1.2  sigma 


resistivity  (ohm-m) 


Figure  20.  Minimum  structure  inversion  of  USGS  data  with  5%  random  noise  added. 
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5.4  Silver  Fox  Mine 

The  data  from  this  experiment  were  collected  by  APTI  during  their  1995  Silver  Fox  Mine 
campaign  using  the  HIP  AS  ionospheric  heater  some  40  km  away  from  the  experiment  site.  This 
transmitter  is  similar  to  that  being  constructed  in  the  HAARP  facility.  We  first  analyzed  the 
amplitude  spectra  from  several  sites  and  for  different  frequencies  of  transmitted  signal.  Figure  21 
shows  the  amplitude  spectra  at  site  4  for  a  transmitter  frequency  of  995  Hz.  While  the  spectra 
have  clearly  defined  spectral  lines  (peaks  in  the  spectra),  these  lines  are  consistent  from  site  to 
site  and  for  the  different  transmitter  frequencies.  Hence,  these  spectral  lines  must  be  related  to 
background  signals  (natural  and  artificial)  and  not  to  the  signals  being  transmitted  by  the  HIPAS 
source.  In  fact,  it  was  difficult  to  find  evidence  of  the  HIPAS  transmitted  signals  except  at  a  few 
sites.  Acting  on  a  suggestion  by  APTI,  we  looked  closely  at  site  5,  and  for  a  transmitted 
frequency  of  2500  Hz  we  were  able  to  identify  spectral  lines  associated  with  the  transmitted 
signal.  The  odd  harmonics  of  this  signal  were  especially  strong. 
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Since  the  transmitted  signals  were  in  general  too  weak  to  even  identify,  we  simply 
proceeded  to  process  the  time  series  using  the  robust  processing  algorithms  developed  in 
conjunction  -with  Dr.  Jim  Larsen  ofNOAA  (Larsen  et  al.,  1996).  At  all  sites  and  for  all 
transmitted  frequencies,  we  found  that  the  TM  mode  data  (electric  field  perpendicular  to  the 
tunnel)  were  noisy  and  had  low  coherencies.  The  TE  mode  data  (electric  field  parallel  to  the 
tunnel),  however,  were  of  much  better  quality.  Unfortunately,  it  is  the  TM  mode  data  that  would 
be  most  sensitive  to  a  resistive  sub-surface  tunnel.  An  example  of  the  results  of  the  robust 
processing  is  shown  in  Figure  22  for  station  15.  In  this  figure,  the  El  direction  is  parallel  to  the 
tunnel  (TE  mode),  and  the  E2  direction  is  perpendicular  to  the  tunnel  (TM  mode). 

The  robust  processed  data  were  then  used  in  a  minimum-structure  inversion.  Both  TM 
and  TE  mode  data  were  included,  but  the  TM  mode  data  were  assigned  larger  variances  to  reflect 
their  lower  coherencies.  The  result  of  the  inversion  is  shown  in  Figure  23  as  a  color-shaded 
cross-section  of  the  subsurface  resistivity.  An  area  of  increased  resistivity  is  indicated  in  this 
figure,  Avith  its  peak  located  around  the  4*  station  from  the  left,  and  at  around  20  meters  in  depth, 
which  is  close  to  the  actual  location  of  the  subsurface  tunnel.  It  should  be  noted  that  the  ability  to 
detect  and  image  this  tunnel  is  not  related  to  the  transmitted  HIP  AS  signals,  which  were  mostly 
too  weak  to  be  recorded,  but  rather  simply  to  the  natural  background  AMT  signals.  If,  however, 
the  transmitted  HIPAS  signals  were  stronger,  they  could  play  a  valuable  role  in  imaging  the 
Silver  Fox  Mine  tunnel. 
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Figure  22.  An  example  of  robust  processing  applied  to  the  data  from  station  15.  The  El 
direction  is  parallel  to  the  tunnel,  and  the  E2  direction  is  perpendicular  to  the  tunnel. 


Figure  23.  Minimum  structure  inversion  result  of  Silver  Fox  Mine  data. 
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6.  Robust  Processing  Algorithm 

AMT  data  present  special  problems,  compared  with  MT  data,  due  to  local  lightning  sources,  low 
signal  levels  around  1-2  kHz,  changes  in  the  transmitter  frequency  between  data  sections  and 
single  site  data  limiting  the  analysis  to  least  squares  biased  by  magnetic  noise. 

The  robust  processing  code  used  to  obtain  transfer  functions  was  originally  developed  for 
MT  data  (Larsen  et  al.,  1 996).  It  was  recently  modified  by  Dr.  Larsen  to  allow  for  the  changing 
frequency  of  the  transmitter  by  developing  a  data  selection  process  for  each  frequency  band. 
Local  lightning  effects  were  excluded  by  examining  the  time  domain  residuals.  It  was  necessary 
to  set  the  outlier  detection  threshold  higher  than  for  MT  data  so  that  only  the  really  local 
lightning  were  excluded.  An  example  of  this  is  shown  in  time  series  collected  in  an  AMT 
experiment  over  a  tunnel  near  Ridgecrest,  California.  In  this  example,  there  is  also  a  local 
induction  transmitter  located  some  200  meters  away  from  the  survey  line.  Lightning  effects  are 
included  in  the  analysis  for  Figure  24  (they  are  not  tagged  as  outliers)  and  excluded  for  Figure  25 
(they  are  tagged  as  outliers  and  removed).  The  apparent  resistivities  and  phases  are  shown  in 
Figure  26  for  the  case  when  the  outliers  are  removed.  The  response  represents  significant 
improvement  over  our  initial  robust  analysis  of  the  data. 

Significant  improvements  in  the  robust  code  are:  (1)  smooth  estimate  of  the  transfer 
function  based  on  a  least  squares  spline  fit  to  the  band  averaged  estimates.  The  smooth  estimate 
is  used  for  estimating  the  residuals  that  are  then  used  to  detect  and  remove  outliers  in  the  time 
domain  and  to  estimate  the  frequency  domain  whitening  and  weights,  (2)  section  weights  based 
on  coherence,  signal/noise  level  and  discrepancy  of  the  transfer  function  from  its  median  value 
for  both  the  entire  frequency  range  and  by  frequency  bands,  and  (3)  smooth  estimates  excluding 
band  estimates  having  low  coherence  or  large  errors.  This  eliminates  the  dead  band  effect  on  the 
smooth  estimates. 

The  robust  method  depends  on  correcting  the  time  series  for  the  filters  using  a  spline 
interpolation  between  the  measurement  filter  response.  This  interpolation  did  not  work  well  for 
the  present  data  due  to  notch  filters  that  cause  the  filter  response  to  change  rapidly  in  phase  with 
frequency.  We  found,  however,  that  the  ratio  of  the  electric  response  over  the  magnetic  response 
is  smoothly  varying  in  frequency  and  can  therefore  be  spline  interpolated.  Correcting  the 
magnetic  series  for  the  ratio  of  the  filters  removes  the  filter  effects  from  the  transfer  function. 
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Figure  24.  Example  of  time  series  containing  lightning  effects  that  were  included  in  the 
determination  of  the  transfer  function.  The  small  periodic  variations  are  due  to  the 
transmitter. 
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Figure  25.  Example  of  time  series  containing  lightning  effects  that  were  excluded  from  the 
analysis. 
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Figure  26.  Transfer  function  for  the  -Zen  element.  Note  that  the  high  and  low  band 
estimates  fit  nicely  together. 

7.  Conclusions 

We  have  reported  on  our  research  from  a  one  year  contract  with  the  Air  Force  Research 
Laboratory  for  advanced  imaging  of  underground  structure  using  HAARP.  Our  research  has 
included  the  following  areas:  1)  adaptation  of  our  minimum  structure  inversion  algorithms  to 
solve  for  simple  parametric  models,  2)  analysis  of  the  sensitivity  and  resolution  of  surface  EM 
data  to  subsurface  tunnels,  3)  analysis  of  EM  data  collected  independently  over  known  tunnels, 
and  4)  modification  of  our  robust  processing  algorithms  to  deal  with  AMT  data  and  controlled 
source  transmitters.  Our  parametric  inversion  algorithms  include  one  to  invert  for  sharp 
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boundaries  of  variable  geometry  and  one  that  inverts  for  a  rectangular  body  in  a  background 
media.  Both  have  proven  successful  for  simple  models.  Our  analysis  of  sensitivity  and  resolution 
of  surface  EM  data  indicates  that  tunnels  need  to  be  at  depth  to  diameter  ratios  of  approximately 
3:1  in  order  to  be  detectable  by  surface  EM  methods.  This  was  borne  out  by  analysis  of  actual 
field  data,  all  of  which  was  of  marginal  quality.  Low  signal  levels  and  noise  resulted  in  noisy  and 
biased  data,  making  accurate  interpretations  difficult.  We  have  modified  our  robust  processing 
algorithms  to  deal  with  AMT  data  and  controlled  sources,  with  a  resulting  improvement  in  the 
responses  in  cases  of  low  signals  and  high  noise  levels. 
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